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Phase structure of monolayer graphene is studied on the basis of a U(l) gauge theory defined on 
the honeycomb lattice. Motivated by the strong coupling expansion of U(l) lattice gauge theory, 
we consider on-site and nearest-neighbor interactions between the fermions. When the on-site 
interaction is dominant, the sublattice symmetry breaking (SLSB) of the honeycomb lattice takes 
place. On the other hand, when the interaction between nearest neighboring sites is relatively 
strong, there appears two different types of spontaneous Kekule distortion (KD1 and KD2), without 
breaking the sublattice symmetry. The phase diagram and phase boundaries separating SLSB, KD1 
and KD2 are obtained from the mean-field free energy of the effective fermion model. A finite gap 
in the spectrum of the electrons can be induced in any of the three phases. 
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I. INTRODUCTION 

The discovery of graphenei, a one-atom thick material 
of carbon atoms, has made a great impact not only on 
condensed matter physics but also on particle physics 2 . 
It gives a realization of massless quasiparticles in a ma- 
terial easy to create and observe, since the valence band 
and the conduction band of the electrons touch at two in- 
dependent "Dirac points" in the Brillouin zone with the 
conical shape- . Due to this "Dirac cone" structure, there 
can be seen several unconventional behaviors character- 
istic to monolayer graphene, such as the high mobility 
of charge carriers and the half-integer quantum Hall ef- 
fect. Since these charged quasiparticles obey the Dirac 
equation around half filling, they are described as mass- 
less Dirac fermions in the (2+l)-dimensional plane, as 
an effective field theory^. Such an effective field descrip- 
tion also has some connections to the high energy physics 
side, such as lattice fermion formulation^, deformation- 
induced gauge fields 6 - - — and the existence of vortex zero 
modes 9 -^. 

The effect of the Coulomb interaction between elec- 
trons is one of the most important problems in graphene 
physics^. Since the Coulomb interaction strength in 
graphene is effectively enhanced by the inverse of the 
Fermi velocity v F ~ c/300 from the ordinary quantum 
electrodynamics (QED), it is beyond the treatment of 
perturbative expansion unless the interaction is screened 
by dielectric substrates such as silicon oxides (Si02). If 
the interaction is sufficiently strong, the electron and hole 
may form an exciton condensate, which may give a finite 
gap in the band structure of graphene. This scenario is 
analogous to the dynamical mass generation of fermions 
in strongly coupled gauge theories such as quantum chro- 
modynamics (QCD), where the spontaneous breaking of 
the chiral symmetry leads to the dynamical mass gap of 
the fermionsi 2 -. In the effective field theory of graphene, 
the chiral symmetry of the fermions corresponds to the 
inversion symmetry between two triangular sublattices 
of the honeycomb lattice. Owing to such an analogy, 



there have been several studies on the "chiral symme- 
try breaking" in monolayer graphene with the techniques 
commonly used in the studies on QCD. Schwinger-Dyson 
equatior>i£~— , 1/N expansio n 16 ' 17 , and the exact renor- 
malization group approach^ have been applied to the ef- 
fective field theory of monolayer graphene. Monte Carlo 
simulations of the effective square lattice model have 
been performed to obtain the critical value of the cou- 
pling constant and the equation of state around the crit- 
ical poin t 19 ' 20 . The author has treated the system as a 
strongly coupled U(l) lattice gauge theory by the strong 
coupling expansion, which is one of the methods to in- 
vestigate the non-perturbative features of the strongly 
coupled gauge theories like QCD^ir— , and has obtained 
the behavior of the (pseudo-)Nambu-Goldstone mode re- 
lated to the chiral symmetry breaking in the low energy 
region . 

In graphene, however, there may be other ordering pat- 
terns than the sublattice (chiral) symmetry breaking that 
may open a finite spectral gap, due to the honeycomb lat- 
tice structur o 25 ' 26 . Kekule distortion, which is described 
by the alternating pattern of the bond strengths like in 
the benzene molecule^ 7 ., is one of those ordering patterns 
without breaking the sublattice symmetry. It can be in- 
duced externally by the effect of some substrates 2 ^ or 
adatoms on the layer—. In the author's previous work, 
it has been found that sufficiently large external Kekule 
distortion may restore the sublattice symmetry which has 
been spontaneously broken in the strong coupling limit 
of the Coulomb interaction^. On the other hand, there 
has been an argument that the Kekule distortion may 
appear spontaneously as a result of the electron-electron 
interaction 9 -^. It is currently a great challenge what or- 
der may appear in the vacuum-suspended graphene due 
to the effectively strong Coulomb interaction. In order to 
treat the ordering patterns characteristic to the honey- 
comb lattice exactly, the analysis of the effective field the- 
ory preserving the honeycomb structure is required^ 2 - - — . 

In this paper, we investigate the competition between 
two phases, the sublattice symmetry broken (SLSB) 



phase and the Kekule distortion (KD) phase, by using 
an effective fermion modei of graphene keeping the origi- 
nai honeycomb lattice structure. This model, motivated 
by the strong coupling expansion of U(l) lattice gauge 
theory defined on the honeycomb lattice, includes the 
on-site interaction and the nearest-neighbor (NN) inter- 
action. If the on-site interaction is dominant, the SLSB of 
the honeycomb lattice takes place. On the other hand, if 
the interaction between nearest neighboring sites is rel- 
atively strong, there appears two types of spontaneous 
Kekule distortion (KD1 and KD2), without breaking the 
sublattice symmetry. By analyzing the mean-field free 
energy of the the effective fermion model, we derive the 
phase diagram and phase boundary separating the three 
phases, SLSB, KD1 and KD2. A finite gap in the spec- 
trum of the electrons can be induced in any of the three 
phases. 

This paper is organized as follows. In Section [Til 
we construct U(l) lattice gauge theory on the honey- 
comb lattice starting from the conventional tight bind- 
ing Hamiltonian coupled with the electromagnetic field 
as compact U(l) link variables. In Section IIII1 we de- 
rive the interaction terms between fermions by using the 
strong coupling expansion techniques of the U(l) lattice 
gauge theory up to the next-to leading order— Two 
characteristic interactions are induced; the on-site inter- 
action which favors SLSB, and the NN interaction which 
favors KD. We take an effective fermion model including 
these two interaction terms. In the next two sections, we 
take the effective fermion model as it is and investigate 
the phase structure by varying the strength of the on- 
site and NN interactions, to study the interplay between 
these different orders. In Section IIV1 we investigate the 
phase structure qualitatively by taking two characteristic 
cases; on-site dominance and NN dominance. Phase di- 
agram with SLSB, KD1 and KD2 phases are also drawn 
qualitatively. In Section |Vl we confirm the the phase di- 
agram in the previous section numerically by minimizing 
the mean-field free energy of the effective fermion model. 
Section I VII is devoted to summary and concluding re- 
marks. 



II. GAUGED HONEYCOMB LATTICE MODEL 

In order to construct the model action of the system 
preserving the original honeycomb lattice structure, we 
start from the conventional tight-binding Hamiltonian 3 , 
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which describes the hopping of an electron between near- 
est neighboring sites with amplitude h = 2.7eV. Here 
a((v) and 6(6^) are the annihilation (creation) opera- 
tors of electrons on the lattice sites in A and B sublat- 
tices respectively. Si = (0, — a), s 2 = f-^ 2 ) , s 3 = 

-^y 2 , 1 ) are the hopping directions, with the lattice 
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FIG. 1: The schematic picture of the Brillouin zone Q and 
the Dirac points K±. When the Kekule distortion pattern 
is induced, the unit lattice in the real space is three times 
enlarged, so that the Brillouin zone is split into three parts, 
Q, and 



spacing a = \si\ = 1.42 A. The triangular sublattices A 
and B are spanned by two lattice vectors Ri = s 2 — Si 
and R2 = S3 — s±. In the momentum space, the Bril- 
louin zone is spanned by reciprocal vectors Kj^, where 
Ki • Rj = 2nSij. By diagonalizing this Hamiltonian in 
the momentum space, the dispersion relation reveals the 
"Dirac cone" structure 



£(K±+k) = fe|4(K± +k)| 



0(k 2 



(2) 



around two independent Dirac points K±, where $(k) = 
2 3 e ik ' Si (see FiglTJ. When the system is half-filled, 
the valence band and the conduction band touches only 
at these points. The Fermi velocity v F — (3/2)ah = 
3.02 x 10~ 3 is considerably smaller than the speed of 
light. This Hamiltonian possesses an inversion symmetry 
between two sublattices A and B, which can be extended 
to the continuous U(1)a symmetry in the low-energy re- 
gion. On-site energy difference between two sublattices, 
m(a'a — b'b), breaks this sublattice symmetry, which cor- 
responds to the mass term rmptp of the Dirac fermions. 

From the Hamiltonian in Eq.(p]), the effective action 
for fermions Sf is derived with the imaginary time (r) 
formulation. Here we perform the temporal scale trans- 
formation r — > t' /v F , so that the Fermi velocity shall be 
rescaled to unity. The temporal direction is discretized 
with the lattice spacing a T /{= v F a T ) equal to the spatial 
lattice spacing a. As a consequence of this discretiza- 
tion, we have a pair of fermion doublers in the temporal 
direction^, which we consider here as the spin (up/down) 
degrees of freedom. 

In this lattice model, the effect of the electromagnetic 
field is implemented by U(l) link variables between spa- 
tially or temporally neighboring sites: 

S F = \ Yl W{x)U Tl {x)a{x + aT') - H.c.j (3) 

r£A;T' 

+- [bHx)U T ,(x)b(x + ar') — H.i 
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where the lattice site x = (r, r') and the link variables 
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, (temporal) (4) 
, (in-plane) (5) 
(out-of-plane) (6) 



This lattice construction is similar to that employed in 
Refl34l. while they differ in the treatment of the U(l) 
gauge field and the imaginary time discretization. 

Using these U(l) link variables, the kinetic term of the 
gauge field 
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S G = ± fd 4 xJ2 idpAv-dvAtf 
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can be rewritten on the honeycomb lattice as 
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where the QED coupling constant g 2 = e 2 /eo = 47to:qed, 
and the constant terms are neglected. The sum is taken 
over the (3+l)-dimensional space. Here, the plaquette 
on the (x, y)-plane, Uhex, is hexagonal-shaped, while the 
other are square-shaped. The plaquettes are defined as 

f7 hex (;r) = Ui(x)U;(x + S! -s 3 )U 2 (x + si - s 3 ) (9) 
xU{(x + s 2 - s 3 )[/ 2 (a; + s 2 - s 3 )[/ 2 (a;) 
U lz {x) = Ui{x)U z {x + Si)Ut[x + az)U* z {x) (10) 
U lT i{x) = Ui{x)U T ,(x + Si)U*(x + aT')U;,{x) (11) 
U ZT ,{x) = U z (x)U T i(x + az)U*(x + a?)U*,{x). (12) 

As a consequence of the temporal scale transformation, 
the spatial part of the gauge field [the first line in 
Eq.flSJ] becomes weakly coupled with the effective cou- 
pling strength g 2 v F , while the temporal part [the second 
line] becomes strongly coupled with the strength g 2 /v F . 
Since the coefficient of the spatial part l/g 2 v F is suffi- 
ciently large, we can apply a saddle point approxima- 
tion to these two terms, yielding a saddle point solution 
Uhcx = Ui Z = 1. In other words, the retardation effect 
(magnetic field) can be neglected due to the discrepancy 
between the speed of light and the speed of fermions (v F ) , 
which is referred to as "instantaneous approximation." 
We can safely set the spatial link variables Ui and U z 
to unity by the gauge transformation, leaving only the 
temporal link variable U T i. The fluctuation of the spa- 
tial link variables around the saddle point, which can be 
considered by weak coupling expansion, is not taken into 



account in this work. As a result, So can be simplified 

as 
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S G = -V3f3 ^2ReU T i(x + Sl )U*,(x) 

r£A;T' i=l 

-3A/3/3 ReU T <[x + az)U;,{x). (13) 

reAUBjT' 

The parameter f3 = v F /g 2 represents the inverse of 
the coupling strength, which is 0.037 in the vacuum- 
suspended graphene. Here, we fix the Fermi velocity v F 
to the physical value observed in the system with Si0 2 
substrate. 



III. STRONG COUPLING EXPANSION 

Here we apply the techniques of the strong coupling 
expansion to the effective action defined in the previ- 
ous section, and derive the effective interaction terms 
between fermions by integrating out the gauge degrees 
of freedom, to construct the effective model which may 
describe the interplay between the sublattice symmetry 
breaking and the Kekule distortion. With the effective 
action S — Sf + So, the partition function of the system 
is given by path integral 

Z = J [d X *dx][dU T ,] exp {~S F [ X \ x; U T ,] - S G [U T ,}} , 

(14) 

where x = a, 6. Since the gauge term Sq is proportional 
to the small parameter j3, we can expand this equation 
around j3 = (strong coupling limit) and perform the 
path integral order by order: 



z = J2 z(n) > ( 15 ) 

n=0 

= ( [dx U x ]\dU T ,}e- s *tMl, (16) 



Here we take the terms up to 0(/3 1 ). Since the integrand 
can be written as a polynomial of U T > and U*, , integration 
by the link variables can be performed analytically. As 
a result of the link integration, two kinds of interaction 
terms are derived: the on-site interaction in the lead- 
ing order [O(/3 )], and the nearest neighbor interaction 
in the next-to leading order [0{j3 1 )]. In order to convert 
these four-Fermi terms into fermion bilinear, we apply 
the Stratonovich-Hubbard transformation by introduc- 
ing two kinds of bosonic auxiliary fields, corresponding 
to the amplitude of sublattice symmetry breaking and 
the spontaneous Kekule distortion respectively. By inte- 
grating out the fermionic degrees of freedom, we derive 
the effective potential of the system as a function of these 
order parameters. 
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Z(r,r'+ai') z f (r,r'+af') a(r,r'+ai') b\r +s J ,r'+ar l ) 




x\r,r') Z(r,r') a\r,'r) b(r + s p r') 

(a) (b) 

FIG. 2: Schematic pictures of the link integration in (a) the 
leading order (LO) and (b) the next-to LO (NLO) in the 
strong coupling expansion, (x = a, b) 



where n x (x) = X i x )x( x ) denotes the local charge den- 
sity at the site x = (r, t'). 

Here we apply Stratonovich-Hubbard transformation 
by introducing the bosonic auxiliary field er, which cor- 
responds to the charge density difference between A and 
B sublattices, (n a — nt). By mean-field approximation 
over er, the first line in Eq. (|20[) is converted into fermion 
bilinears as 
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A. Leading order: O(/3 ) 

In the leading order (LO), the gauge term Sc does not 
contribute, 



Z (a) = [d X U X }[dU T ,}e- s - , 



(17) 



so that the link variables come only from the temporal 
hopping terms of the fermions. On each lattice site x = 
(r,r'), the contribution to the link integration is 



dU T > exp 
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(18) 
(19) 



where x' = x( x + ar'). The schematic picture of the link 
integration in the LO is shown in Figl^a). By the link 
integration, on-site four-Fermi interaction term is gener- 
ated, which corresponds to the on-site repulsion between 
opposite spins (Hubbard term). In order to control the 
strength of this interaction for later purpose, we intro- 
duce an overall coefficient z(> 0). Thus, the effective 
action can be written in terms of fermionic fields a and 
b as 



o(0) z 

F ~~4 



Y, n a(x)n a (x + or') + Y^ n b( x ) n b(x + ar' 



r£A;T' 



reB;r' 



- Y Y [aHx)b(x + s i ) + b^(x + s i )a(x)],{20) 



Thus, we can integrate out all the fermionic degrees of 
freedom, to obtain the effective potential of this system 
at LO per one pair of A and B sites: 



1 



N T >V 



lnZ<°> 



(22) 
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where V is the number of A (B) sites in the system. J kgf2 
is the integration over the Brillouin zone f2, with normal- 
ization y /ken < ^ 2 ' s ' 1 = 1- clN t i is the temporal lattice 
size, corresponding to the inverse temperature. In this 
work we consider the zero-temperature and infinite vol- 
ume limit, so that N T > and V are set to infinity. The first 
term in Eq. (|2"2"|) represents the tree level of <r, while the 
second logarithmic term comes from the one-loop effect 
of the fermion. 



B. Next-to leading order: O^ 1 ) 

Next, we consider the next-to LO (NLO) terms in the 
strong coupling expansion, At O^S 1 ), one plaque- 

tte from So contributes to the link integration. Sq (with 
instantaneous approximation) contains two kinds of pla- 
quettes, Ui T i and U ZT >, but U ZT i does not contribute to 
the link integration, because the link in the z-direction 
cannot be canceled by the fermion hopping terms. On 
the other hand, Ui T >(x) = U T >{x + Si)U*,(x) contributes 
to the link integration, combined with two fermion hop- 
ping terms: 



reAir' i=l 
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^ T ,( x ) rf ^ T ,( x + Si ) e -M« t W^Kx)o(x+ a r')+6 t (»+Si)t/ T ,(x+s i )6( a =+s i +ar')-H.c.] x ^M-U^X + Bi)U^,{x) 



dU T >(x) 



U*,(x) — -^a)a + -a^U*, 2 (x)a - 



—a) a U*,{x)a 



x / dU T '(x + Si) 



U T /(x + a T ,) - -tfU T ,(x + 



a T ,fb' + ifc'tft - \tfU T ,(x + aT>)b'b'% 



a'(x)a(x + ar')w{x + s, + ar')b(x + Sj), 



(23) 



(24) 



as shown in FigJ^b). (In Eq. (j2"3"|) . we denote a = 
a(x), a' = a(x + ar'), = b(x + Sj) and b' = b{x + 
Si + ar').) Thus, the effective action in the NLO can be 
written in terms of fcrmions as 



5' 



(i) 



-*E E 

r£A;r' i=l 



(25) 



^a f (a;)6(x + s, t )tf (x + s, + ar')a(x + a?) + H.c. 

Hereafter, we use the rescaled parameter £ = \/3/3/8 in- 
stead of /3 as the strength of such a nearest-neighbor in- 
teraction. 

In order to convert this interaction term into 
fermion bilinears, we apply the "extended" Stratonovich- 
Hubbard transformation, 



aAB = const, x / d\d\*e- a ^ 2 - XA - x ' B \ (26) 



with the positive constant a and the complex auxiliary 
field A. With the auxiliary field A,(r, r') corresponding 
to the fermion bilinear (<r(x)b(x + Sj)), we obtain the 
NLO effective action 



=2£J2 E 0M*)| 2 - (Ai^at^ft^ + sO+Hx.)] 

reA-r' i=l 

(27) 

Thus, S 1 ^ modifies the hopping of the fermions in the 
spatial direction through the auxiliary field A^. 

Here, we take the ansatz that Aj should be split into 
the spatially uniform part and the spatially varying part 
with the Kekule distortion pattern: 



e i(K + -s 3 +G-r) . ( < K s C; 



Xj(r,T') = X a +X A e 2 ^ 3 

(28) 

where G = K + — K , and A CT and Aa are real values. 
The first term renormalizes the Fermi velocity v F uni- 
formly, with the factor Z v = 1 — £A CT /3. We show later 
that (A<j) < 0, so that the Fermi velocity becomes faster 
at finite j3 (or £) than that in the strong coupling limit 
(/3 = 0). The second term corresponds to the sponta- 
neous Kekule distortion, with the amplitude A = 3£Aa- 




FIG. 3: Schematic picture of the Kekule distortion pattern. 
Thick lines and thin lines represent the strong hopping and 
the weak hopping, respectively, (a) Distortion pattern for 
Aa > 0. (b) Distortion pattern for Aa < 0. 



The Kekule distortion is characterized by the pattern 
of alternating bond strengths, as shown in FigOJ and 
induces a spectral gap without breaking the sublattice 
(chiral) symmetry-, with the modified dispersion rela- 
tion E(K± + k) = Vlkp + fAp + 0(k 4 ). Since its unit 
lattice is three times as large as that of the ordinary hon- 
eycomb lattice in the real space, the Brillouin zone fi is 
split into three hexagonal parts: and Q±, surrounding 
T-point and the Dirac points K± respectively, as shown 
in FigfTJ As a result, the effective action up to the NLO 
can be written with the order parameters a, X a and Aa 
as 



4° ) +4 1) = E [^ 2 + 6^ + 2A A ) 



r£A;r' 



(29) 



E * f ( k < T ') 



-{z/2)ah (2/3)$t( k ) \ - , 
(2/3)$(k) (z/2)al 3 ) 1 ' h 



where the 3 x 3 matrix $(k) = $o(k) - 3£<li(k), with 
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A' A $(K_+k) A CT $(K++k) A A $(k) 
v A' A $(K++k) A A $(k) A ff $(K_+k) 
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and I3 is a 3 x 3 unit matrix. Here we denote A' A = 
AAe~ 27 ™/ 3 for simplicity. The fermionic field ^> is de- 
fined by *(k,r') = [a(k, r'), a(K+ + k,r'),a(K_ + 

k, r'), 6(k, r'), &(K+ + k, r'), 6(K_+k, r')] T . By integrat- 
ing out the fermion field ^, we obtain the effective po- 
tential 



F 



(0+1) 



cff {a, \ a , Aa) = 2^ 



V 



d 2 kln det 



kef2 



2 + 6£(A 2 

(?)V 



2A|) 



(31) 



3) $ f (k)$(k) 



The third term (fermion loop effect) is modified from that 
in Eq. (|22]) by the spontaneous Kekule distortion Aa- 



IV. QUALITATIVE PROPERTIES 

So far we have reconstructed the effective fermion 
model with two interaction terms, the on-site interaction 
and the NN interaction, obtained by the strong coupling 
expansion of the U(l) lattice model, and derived the ef- 
fective potential of the system: F^ +1 \a, X^, Aa)- Here- 
after, we vary the strengths of these interaction terms (z 
and £ respectively) to arbitrary values, to observe the in- 
terplay between the sublattice symmetry broken (SLSB) 
phase and the Kekule distortion (KD) phase. First we 
investigate the qualitative properties of possible phases 
in the system by taking the characteristic limits of z and 
£: the SLSB phase in the limit £ ~ 0, and the sponta- 
neous KD phase in the limit z — 0. Then, we consider 
the competition between these two phases by approxi- 
mating the effective potential in the region where both 
z and £ are considerably small, and estimate the phase 
structure of the system qualitatively. As a result, we find 
that the appearance of the SLSB phase or the KD phase 
is related to the dominance of the on-site term or the NN 
term respectively, and that the KD phase is split into two 
phases (KD1 and KD2), flipping the sign of Aa- 



A. Sublattice symmetry broken phase: £ ~ 

First we consider the limit £ — 0, where only the on-site 
interaction is concerned. In this limit, the effective poten- 
tial in Eq. (j3"T1) reduces to the simpler one in Eq. ([2"2"j) . The 
first term (tree level of a) becomes dominant as \a\ — > 00, 
while the second term (fermion one-loop effect) domi- 
nates when \a\ — > 0. Due to the logarithmic singularity 
of d 2 F^ /da 2 around a = from the one- loop term, 

F^\a) has a minimum at finite a for any value of z > 0. 
The potential minimum gives the expectation value of 
the charge density imbalance between two sublattices, 
(a) = (n a — rib), which serves as the order parameter 
of the spontaneous sublattice (chiral) symmetry break- 
ing. In the 4-component Dirac fermion representation, 
it corresponds to the chiral condensate Therefore, 



0.000 
-0.002 




FIG. 4: The behavior of the effective potential of the system 
in the strong coupling limit, F c ^' , as a function of exciton (chi- 
ral) condensate a at z = 1. Honeycomb: F^(a) in Eq. (|22[) 
with the exact dispersion relation "^(k) = ~Y^ i e^' Si . Linear: 
F^ (a) in Eq. ()22p with the approximate dispersion relation 
$(K±+k) = ^e- 2wi/3 a(±k x + ik y ). Square: The effective po- 
tential obtained from the square lattice formulation [Eo. (l32l) ]. 



the sublattice symmetry of the system is spontaneously 
broken in the limit £ = 0. 

The behavior of the effective potential in Eq. ([22| at 
z = 1 is shown in FigfJ] as the curve with the label 
"Honeycomb." In FigfJ] F^\a) obtained from two other 
formulations are displayed: "Linear" is the effective po- 
tential calculated with the Dirac cone approximation 
$(K± + k) = %e- 2 ™^a{±k x + ik y ), and "Square" is 
the one obtained from the square lattice formulation 24 , 



F 



ctf ( CT ) = y 



1 



(»* 



dk 2 In 



ke[- 



(f) 1 



sin 2 kn 



(32) 

All of them qualitatively have the same structure because 
they have the Dirac cone structure in common around the 
Dirac points, but quantitative behaviors are different due 
to the deviation from the Dirac cone structure at large 
momentum. 

As seen from Eq. (|2~Tj) . finite a induces an on-site energy 
difference between two sublattice in the sense of mean- 
field, yielding a finite spectral gap E(K±) = v F za/2a, 
which corresponds to the dynamically generated mass of 
the fermion. When z takes the physical value z = 1 , the 
expectation value is a — 0.343, which gives the dynam- 
ical gap 0.72eV. By taking the momentum integration 
around k = 0, we have an approximate relation 



OF. 



(0) 



off 



da Z + 47r(7 V2 



In 



(33) 



for sufficiently small z and a. Thus, z-dependence of (a) 
around z can be approximated as 



(2/z)exp(-2/z), 



(34) 



so that (a) reaches toward as z 0. 

Next, we observe the behavior of the effective potential 
in the vicinity of £ = 0. In order to take the effective 
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potential up to O^ 1 ), we have to expand the fermion 
determinant by £. By using the formula 

lndetL4 + £B] = In dct A + Tr^- 1 ^] + 0(£ 2 ), (35) 

the third term in Eq. (|3"Tj) is approximated up to O^ 1 ) as 



(?)" 



In det 



lndet 



-3£Tr 



Since the matrices 



(*£-3e*I)($ -3£$i) 



3) *6*o 



3 ) mo 



(36) 



(za/2) 2 + (2/3) 2 Sp 



and $0 are 



diagonal, only the diagonal part of $1, which can be 
written as A CT $o, contributes to the trace in the sec- 
ond term in Eq. (|36l) . Thus, we have the ^-expansion of 

?(0) 



the effective potential as F^ +1 \a, A CT , Aa) 
F e W(a,A CT ,A A ) + 0(C 2 ), where 

F$(a, A CT ,A A ) =6^(A 2 +2A A ) 



(37) 



6gA 

1/ 



d 2 k 



(2/3) 2 |$(k)|^ 



ken 



(z ( x/2) 2 + (2/3) 2 |$(k)| 1 



Taking the potential minimum by the NN-auxiliary 
fields A CT and Aa, we obtain their expectation values up 
to the LO, 



Aa = O + O^ 1 ) 



A„ 



2V 



(2/3) 2 |$(k)| 5 



ken 



(z < r/2)> + (2/3)2|*(k)|s 



(38) 

o(e). 



Since Aa does not contribute to the fermion one-loop 
term up to O^ 1 ), the Kekule distortion does not ap- 
pear around the limit £ = 0. On the other hand, A CT 
acquires a negative expectation value, so that the renor- 
malization factor of the Fermi velocity, Z v = 1 — 3£A ff , 
becomes larger than unity. By substituting these rela- 
tion to F^g (a, A ff , A A ), the NLO effective potential can 
be rewritten as a function only of a: 



d 2 k 



(2/3) 2 |$(k)|; 



ken 



(W2) 2 + (2/3) 2 |a>(k)|^ 



. (39) 

Since this term monotonically increases as a function of 
<r, it reduces the expectation value of a (that is, the po- 
sition of the potential minimum) . At the physical value 
z = 1, the expectation value of a is given up to NLO as 



a(i) = 0.342- 1.73£ + 0(£ 2 ), 



(40) 



and Act = —0.471 + 0(£). Therefore, the system reveals 
the SLSB phase in the vicinity of £ = 0, and the am- 
plitude of SLSB (charge density imbalance between two 
sublattices), a, decreases as a function of £. 



B. Kekule distortion phase: z — 

In order to investigate the qualitative properties of the 
Kekule distortion (KD) phase, we take the limit z — 0, 
where the system does not contain the on-site interaction 
so that it may not reveal the SLSB phase. Here, the 
effective potential reads 



^ +1) (A.,A A ) 



= 6£(A 2 + 2Ai) - - 
= 6£(A 2 + 2A A ) 



<T kin dct 



In 



2 - Zl 
3 



- 2(2£A A ) 



ken 

r 

ken 
3 



!*(k) 



d 2 k- 



(41) 
(42) 



$(k)$(K 4 



X(2£A A ) 2 e 2 ™ /3 [$ 3 (k) + $ 3 (K H 



-k)$(K_ +k) 
k) + $ 3 (K_+k)] 



First we analyze whether Aa takes a finite expec- 
tation value or not. Since one can easily see 

9F c ^ +1 ^/9Aa|a a =o = 0, Aa = is either a local max- 
imum or minimum of the effective potential. In order to 
consider the behavior around Aa, we have to check the 
sign of the second derivative 



d 2 F, 



(0+1) 



A A =0 



= 2A£ + — 
* V 



d 2 k- 



ken 



^ 2 [$ 3 (k) 



3> 3 (K+ +k) + $ 3 (K_ 



(43) 
k)] 



(|Z„) 3 e- 2 W3$(k)$(K H 



k)$(K"_ +k) 



Since the denominator of the integrand becomes zero only 
at k = 0, the region around this point is dominant in 
the loop integration. Taking the leading order in |k| in 
the numerator and the denominator, the loop integral 
becomes 



2 
V 



144£^ 
ken 2Z v \k\ 2 - 



0(|k| 3 ) 



(44) 



which has a negative logarithmic divergence. Due to this 
logarithmic divergence in the momentum integration, the 
sign of the second derivative becomes negative at Aa = 
0, so that Aa = is a local maximum of the effective 
potential. Therefore, for any value of £(> 0) (or /?), Aa 
takes a finite expectation value. 

Next, we consider the behavior of the potential mini- 
mum (A CT , Aa). Due to the logarithmic divergence of the 
momentum integration when the order parameters sat- 
isfy the relation § - 2£A CT - 2{/2^X A = [see Eq.(@2])], 
the effective potential has a non-analyticity on this curve, 
separating the (A CT , AA)-plane into two regions. In each 
region there is a local minimum of F c g, and it depends on 
the value of £ which minimum is taken. When £ crosses 
over a certain value , one local potential minimum may 
dominate over the other one, causing a sudden jump of 
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(Aa). Therefore, the KD phase is split into two regions 
at the line £ = where the system reveals the first 
order phase transition. Here we refer to these two phases 
as KD1 for !; < £k and KD2 for £, > £,k, respectively. 

Finally, we consider the behavior of (Aa) in the limits 
£ ~ and £ — > oo. In the limit £ ~ 0, we neglect the A CT - 
dependence because it depends on the loop integration 
only via Z v = 1 — 3£A CT , which becomes unity at £ = 0. 
By performing the momentum integration around k = 0, 
we have the approximate relation 

F^ +1) (A A ) ~ 12£A A + A(2£\ A ) 2 m(2£A A ) 2 , (45) 

where A is a positive constant related to the area of the 
momentum integration. By taking the potential mini- 
mum, we can estimate the order of (Aa) to be 

(A A )^r 1 cx P (-r 1 ). (46) 

Therefore, (Aa) reaches toward zero as £ — > 0. On the 
other hand, in the limit £ — > oo, all the terms in | ■ • • | 
in Eg. (14"2"T) becomes proportional to £ 3 . Thus, the in- 
dependence in the logarithm can be factored out, so that 
only the first term (tree level of A CT and Aa) becomes dom- 
inant in this limit. Therefore, both the expectation val- 
ues of A a and A CT reach toward zero in the limit £ — > oo. 

C. Competition between SLSB and KD phases 

Let us now investigate the competition between two 
phases, SLSB and KD, and observe what kind of phase 
transition may occur between these two phases. In order 
to treat the logarithmic singularity of the loop integral, 
we take into account the momentum space only around 
the Dirac points. Here we consider the region where the 
interaction strengths z and £ are sufficiently small, to 
simplify the discussion. Since £Aa reaches toward zero 
as £ — > 0, we assume that the terms of 0(£A A |k|) and 
the smaller ones are negligible, which gives a simplified 
form as follows: 



^ diag {(? 



h + [-) $T(k)$(k) 



AZi 



2 

(^)%|^k| 2 + 36(£AA) 2 }. (47) 

This simplification is valid as long as the logarithmic sin- 
gularity is dominant, that is, a and £Aa are in vicin- 
ity of 0. Since the first element of this matrix does not 
contribute to the logarithmic singularity, we can further 
simplify this model by neglecting the first element (con- 
tribution from the Brillouin zone 17, which does not cover 
the Dirac points K±). Thus, we obtain the effective po- 
tential 



^ +\z v k\ 2 + 



36(£A A ) 2 , 



^ +1) (M ff ,A A )^fa 2 



6£(A 2 



2A A ) 



(48) 



2 
V 



d 2 kln 



ken 



(f) +»M 



\z„k\ 



The properties of the effective potential in Eq. (|48|) can 
be observed rather easily than the exact one. If we define 
a new field <j> by 4> 2 = a 2 + (144£ 2 /.z 2 )A 2 1 , the effective 
potential is rewritten as 



, Ao-, Aa) 



f^ 2 +12dl-f) . • (4!.. 



2 

* J ken 



d 2 kln 



If 1 — 6£/z = (z = 6£), the effective potential is given 
as a function of <j> and X ai and does not depend on Aa 
explicitly. Therefore, a and Aa can take arbitrary expec- 
tation values satisfying 



144£ 5 



(50) 



If 1 - 6£/z > (z > 6£), the second term in Eq.(@9j be- 
haves as a symmetry breaking term between a and Aa- 
Since the fermion loop integral does not explicitly depend 
on Aa, the effective potential monotonically increases as 
a function of Aa, yielding (Aa) = 0. On the other hand, 
4> contributes to the loop integral, so that ((f) ^ 0. There- 
fore, we have (a) =/= 0, that is, the sublattice symmetry 
is spontaneously broken. 

If 1 — 6£/z < (z < 6£), the coefficient of the sec- 
ond term in Eg. (14T)1) becomes negative, leading to the 
unphysical result (Aa) = oo. Here we rewrite the effec- 
tive potential as a function of <p,a and <p a , so that all the 
coefficients of these variables would be positive: 



F cS ((f>,a, X a ) = — > 
2 

~V 



z I z 
2 I ~ 6£ 



a 2 +6£A 2 (51) 



d 2 kln 



keo 



\ZM 



In this form, the effective potential monotonically rises 
as a function of a, so that we have (a) = 0, (<j>) ^ 
and (Aa) ^ 0. In other words, there appears a Kekule 
distortion pattern spontaneously. 

Therefore, when crossing the line z = 6£, there is a 
first-order phase transition between the sublattice (chi- 
ral) symmetry broken (SLSB) phase and the spontaneous 
Kekule distortion (KD) phase. Moreover, as shown in 
the previous subsection, the expectation value of Aa re- 
veals the non-analyticity at a certain value £k in the KD 
phase. Since the effective potential does not depend on z 
in the KD phase, the value of is independent of z as 
long as the point (z,£k) is in the KD phase. Since the 
KD1 and the KD2 phases correspond to different poten- 
tial minima respectively, the critical line between SLSB 
and KD1 and that between SLSB and KD2 are discon- 
tinuous. From the qualitative discussions above, we can 
map a schematic phase diagram of the system, as shown 
in FigEl 
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(c) 
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FIG. 5: A hypothetical phase diagram of the effective model 
of monolayer graphene, from the qualitative investigation of 
the effective potential i^2 + • When the on-site interaction 
(z) is dominant over the nearest-neighbor interaction (£), the 
system spontaneously shows the sublattice symmetry broken 
(SLSB) phase; otherwise the system shows the Kekule dis- 
tortion (KD) phase. The KD phase is separated by the line 
£ = £x into two regions, which correspond to the two local 
minima of the effective potential respectively. All the critical 
line show the first-order transition behavior. The dashed lines 
(a), (b) and (c) correspond to the results from the numerical 
calculation: (a) corresponds to z = 1 in Fig[7] (b) to z = 40 
in FigOO and (c) to z = 50 in FiggO 



V. NUMERICAL RESULTS 

Now we confirm the qualitative results above by min- 
imizing the exact effective potential in Eg. pip numeri- 
cally. In the limit 2 = 0, the effective potential becomes 
independent of a, so that we only derive the expectation 
value of Aa as a function of £, as shown in FigJSJ It 
can be clearly seen that Aa(£) shows a non-analyticity 
at £,k = 9.97, which separates the KD phase into two 
regions. In the KD1 region (£ < £k), Aa obtains a pos- 
itive expectation value, which corresponds to the lattice 
distortion pattern shown in Figf5Ja). As qualitatively 
estimated, (Aa) starts from zero at £ = and monoton- 
ically increases for small £ . It has a peak at f ~ 0.7 and 
eventually decreases until £ = On the other hand, in 
the KD2 region (£ > Aa takes a negative expecta- 
tion value, corresponding to the pattern in Figj3^b). It 
then reaches toward zero as a function of £. which agrees 
with the analytical result that (Aa) — > as £ — >• oo. 

Next, we fix the on-site interaction strength z to fi- 
nite values. At the value z = 1, which corresponds to 
the strong coupling expansion of the gauged model, the 
expectation values of the sublattice symmetry breaking 
amplitude a and the spontaneous Kekule distortion Aa 
vary as function of £, as shown in FigJT] The order pa- 
rameters reveal non-analyticity at two points, £ = 0.20 
and £ = In the region £ < 0.20, a is finite and 

monotonically decreases, while A a is zero. As can be 
clearly seen from the analytic observation, this region 
corresponds to the sublattice symmetry broken (SLSB) 




10 12 14 



FIG. 6: The behavior of the spontaneous Kekule distortion 
(KD) Aa as a function of the nearest-neighbor interaction 
strength f, in the absence of the on-site interaction (z = 0). 
The inset shows the behavior of Aa in vicinity of £ = 0. Aa 
shows a non-analyticity at £k = 10.67, where the KD phase is 
separated into two regions, KD1 (Aa > 0) and KD2 (Aa < 0). 
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FIG. 7: Expectation values of the spontaneous sublattice (chi- 
ral) symmetry breaking a and the Kekule distortion ampli- 
tude Aa, calculated at the physical value 2 = 1. The inset 
shows the behavior of a and Aa in vicinity of £ = 0. There 
is a first order phase transition from the SLSB phase into the 
KD1 phase at £ = 0.20, and that from KD1 into KD2 at 
£k = 10.67. Such a behavior corresponds to the line (a) in 
the phase diagram in Fig[5] 



phase in the hypothetical phase diagram in FigJS] In 
the region 0.20 < £ < the expectation value of a 
vanishes, while Aa acquires a positive expectation value, 
which corresponds to the KD1 phase. For large value 
of £, Aa shows a non-analyticity at £ = obtains a 
negative expectation value, and reaches toward zero just 
as seen in the z = limit. This region corresponds to 
the KD2 phase. Therefore, we can conclude that the axis 
z = 1 corresponds to line (a) in FiglU that is, the system 
turns from SLSB into KD1 at a certain critical value of £ 
(here £ = 0.20) and turns from KD1 into KD2 at £ = £ K . 

When z — 40, there appears only one phase boundary, 
as shown in FigjSJ The phase transition occurs at £ = £k , 
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FIG. 8: Expectation values of the spontaneous sublattice 
(chiral) symmetry breaking a and the Kekule distortion am- 
plitude Aa, calculated at z = 40. There is a first order 
phase transition from the SLSB phase into the KD2 phase 
at £jf = 10.67, and the KD1 phase does not appear. Such a 
behavior corresponds to the line (b) in the phase diagram in 

FigG2 




FIG. 10: The phase boundary between the sublattice sym- 
metry broken (SLSB) phase and the Kekule distortion (KD) 
phase, as a result of the numerical calculation. As can be seen 
from Fig[6] the KD phase is split into two phases by the line 
£ = 10.67(= KD1 (Aa > 0) and KD2 (Aa < 0). This 
phase diagram agrees with the qualitative estimation obtained 
in FigEl 
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FIG. 9: Expectation values of the spontaneous sublattice 
(chiral) symmetry breaking a and the Kekule distortion am- 
plitude Aa, calculated at z = 50. There is a first order 
phase transition from the SLSB phase into the KD2 phase at 
£ = 11. 45 (7^ £k ), and the KD1 phase does not appear. Such 
a behavior corresponds to the line (c) in the phase diagram 
in FigE 



from the SLSB phase (a 7^ 0) into KD2 phase (Aa < 
0), and the KD1 phase does not appear. This behavior 
corresponds to the line (b) in FigJSJ 

At the value z = 50, there appears only two phases 
as observed at z = 40, but here the critical value of £ 
is shifted from = 10.67, as shown in FigJHl This 
behavior corresponds to the axis (c) in FigJSJ 

Finally, we show in Fig llOl the phase boundary between 
the SLSB and the KD phases, which agrees with the qual- 
itative phase diagram obtained in FigJSJ In general, the 
system turns into the KD phase when the nearest neigh- 
bor interaction (£) becomes dominant over the on-site in- 
teraction (z), that is, the Coulomb interaction strength 
becomes smaller. Since the effective potential becomes 



independent of z in the KD region, there is a first or- 
der phase transition between KD1 and KD2 at the line 
£ = £k, as seen in the z = limit. 



VI. CONCLUSIONS AND OUTLOOK 

In this work, we have investigated the possible phase 
structure of monolayer graphene with the on-site and 
the nearest neighbor (NN) interactions between fermions. 
First, the effective action of the system is constructed in- 
cluding the electromagnetic field as U(l) link variables, 
and the interaction terms between fermions are derived 
by applying the techniques of the strong coupling expan- 
sion of the lattice gauge theory. Thus we have obtained 
two kinds of effective interaction terms: the on-site inter- 
action which may contribute to the sublattice symmetry 
breaking (SLSB), and the NN interaction which may lead 
to the Kekule distortion (KD). Using these two interac- 
tion terms, we have reconstructed an effective model of 
graphene with arbitrary interaction strengths z and £ re- 
spectively, to investigate the interplay between the SLSB 
and the KD. We have observed the behavior of the order 
parameters with this effective model, by the mean field 
approximation over the effective potential. 

Focusing on the logarithmic singularity of the effective 
potential, we have qualitatively obtained the phase di- 
agram shown in FigJSJ When the on-site interaction is 
dominant, the sublattice (chiral) symmetry of the sys- 
tem is spontaneously broken, leading to the dynamical 
mass term of the fermions. On the other hand, when the 
nearest-neighbor interaction is sufficiently large, the hop- 
ping parameters in the lattice get renormalized with the 
Kekule distortion pattern. In this case the fermions still 
obtain a dynamical spectral gap, without breaking the 
sublattice (chiral) symmetry. Moreover, this KD phase is 
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split into two regions KD1 (Aa > 0) and KD2 (Aa < 0), 
corresponding to two different minima of the effective po- 
tential. Such a splitting line £ = is numerically seen 
by taking the limit where the on-site interaction is omit- 
ted (z — 0). For instance, when the on-site interaction 
strength z = 1, which corresponds to the strong coupling 
expansion of the Coulomb interaction, the system reveals 
the SLSB phase in the strong coupling limit. The system 
turns into the KD1 phase at £ = 0.20 {0 = 0.92), and 
into the KD2 phase at £ = = 10.67. Since z = 1 
and (3 = 0.037 (£ = 0.008) in the vacuum-suspended 
monolayer graphene, we expect a gapped phase with 
SLSB, while the system may reveal the KD phases if the 
Coulomb interaction is suppressed (/3 is increased) by the 
screening effect by substrates or the renormalization of 
the Fermi velocity Up^ 6 -. The KD1 phase does not appear 
at sufficiently large z, as seen at z = 40 and 50 in this 
work. It has been verified both qualitatively and numer- 
ically that all the phase transitions in the phase diagram 
obtained in this work are first order phase transitions, 
that is, the order parameters reveal non-analyticity when 
crossing the phase boundaries. 

There are still several open questions to be solved 
within the framework of this study. Since the spin de- 
grees of freedom are absorbed in the fermion doubling, 
which is the artifact of the lattice discretization, spin- 
related ordering, such as the spin density wave (SDW) 
phas o 37 ' 38 and the "spin-Kekule" phased, cannot be 
identified out of the SLSB and KD phases in this work. 



Some other lattice discretization scheme that exactly 
treats the spin degrees of freedom is needed. The effect 
beyond the NLO is also an interesting issue. The next- 
to-NLO [0(/3 2 )] term, which includes the four-Fermi in- 
teraction between second nearest neighboring sites, can 
spontaneously generate an effective magnetic flux in the 
honeycomb plaquette, leading to the so-called "quantum 
anomalous Hall (QAH)" state 26 . For example, Refl33l 
treats several types of instabilities by the exact renor- 
malization group method on the honeycomb lattice, and 
shows that only four instabilities, SLSB, KD, SDW and 
QAH, may occur by the effect of the Coulomb interac- 
tion, but the competition among these orders is left for 
further investigation. Extension of the lattice strong cou- 
pling expansion method to the bilayer graphene system 
is also required, since a gapped phase has recently been 
observed in bilayer graphene experimentall y 39 ' 40 . Quite 
a rich phase diagram is expected both in monolayer and 
bilayer graphene systems. 
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